Entanglement and parametric resonance in driven quantum systems 
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We study the relationship between entanglement and parametric resonance in a system of two 
coupled time-dependent oscillators. As a measure of bipartite entanglement, we calculate the linear 
entropy for the reduced density operator, from which we study the entanglement dynamics. In 
particular, we find that the bipartite entanglement increases in time up to a maximal mixing scenario, 
when the set of auxiliary dynamical parameters are under parametric resonance. Moreover, we 
obtain a closed relationship between the correlations in the ground state, the localisation of the 
Wigner function in phase space, and the localisation of the wave function of the total system. 
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I. INTRODUCTION 

A quantum system composed of two or more entangled 
subsystems has the intriguing property that although the 
state of the total system can be well defined, it is impos- 
sible to identify individual properties for each one of its 
parts. Understanding the properties displayed by the en- 
tanglement of physical systems is one of the fundamental 
purposes of quantum information theory [l|, Q • In this 
paper we focus on the relationship between entanglement 
and parametric resonance in a composite quantum sys- 
tem. 

In classical dynamics, parametric resonance can occur 
when an appropriate parameter of a system is varied pe- 
riodically in time [||. Under parametric resonance, stable 
points of the undriven system become unstable for spe- 
cific values of the period of the parameter variation. On 
the other hand, an unstable point may become a stable 
one a phenomenon known as parametric stabilization. 
In the quantum regime, it is well known that a periodi- 
cally driven system follows a unitary dynamics. The dy- 
namics is not closed, however, since the external driving 
changes the physical parameters of the system, altering 
its total energy. The usual technique for treating such 
systems is the Floquet Formalism 0, H[ . Perelomov and 
Popov @, and most recently Weigert 0, have studied 
quantum manifestations of classical parametric resonance 
in noncomposite quantum systems in the context of Flo- 
quet theory for the Schrodinger equation. They found, 
for example, that for parameters in the stable regions, a 
discrete spectrum of quasienergies exist. In contrast, in 
the unstable zones, there is a continuous spectrum. 
The question that naturally arises is: what are the quan- 
tum signatures of classical parametric resonance? There 
is no classical behavior analogous to entanglement, it is 
unique to quantum mechanics. In this paper we show 
that this intrinsically quantum property shows strong sig- 
natures of classical parametric resonance in a system of 



two driven coupled harmonic oscillators. Furthermore, 
we find that the entanglement dynamics is related to the 
dynamical behaviour of the Wigner function and the to- 
tal system's ground-state wave packet. This paper is or- 
ganized as follows: In section [TT] we describe the model. 
In section IIIII we introduce the required basic aspects 
of the Lewis-Riesenfeld canonical method and the single 
parametric oscillator. In section IIVI we find the exact 
solution of the Schrodinger equation and discuss the sta- 
bility properties of the auxiliary dynamical parameters 
and their relation to the phenomenon of parametric res- 
onance. In section [V] we present the calculation of the 
entanglement dynamics and the Wigner function. Fi- 
nally, discussion of the results and conclusive remarks 
are presented in sections I VII and I VIII respectively. 



II. MODEL 

Time-periodic quadratic Hamiltonians arc known to 
provide examples of parametric resonance. This type of 
Hamiltonian has the general form 



I 2 

H(t) = - X>iW$ + ) + 7(*)*i5a- 



(1) 



In this work we are interesting in studying the particular 
model described by the time-periodic Hamiltonian 



! 2 

H(t) = - ]T(p 2 + <^i) + 7(t)xix 2> 



(2) 
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where j(t) = g + Ag cos Qt, to is a time-independent fre- 
quency, g is called the static coupling, and Ag is a frac- 
tion of g. 

Further motivation for studying this type of Hamilto- 
nian comes from the time-dependent Dicke model. In 
1954, Dicke showed, in a now celebrated paper Q, the 
coherent spontaneous emission arising from many atoms 
emitting collectively. The single mode Dicke Hamiltonian 
models the interaction of N atoms with a single mode 



2 



bosonic field via dipolc interactions within an ideal cav- 
ity. The time-dependent generalization of Dicke model is 
described by the Hamiltonian 



H{t) = u>o(t)J z 



-uj(t)^a+^i{J + +J-) 



(3) 



By considering the Holstcin-Primakoff transformation of 
the angular momentum algebra and taking the limit 
j — ► oo in the Hamiltonian Eq. ([3]), one obtains an 
effective Hamiltonian 0, [nj, which corresponds to the 
Hamiltonian of two time-dependent oscillators linearly 
coupled through the interaction strength X(t). In terms 
of the field coordinate x and the atoms coordinate y, this 
Hamiltonian becomes 



H(t) 



W (U (-1 , 2»2\ , w o(*) 1-2 , 2«2\ , 

2A(i) v /< 7 ^o xy- E , (4) 



where w(io) 



w (io) 



w 0; 



and £?o(*) 
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+j'wo(i)- In the next section we describe the Lewis- 
Riesenfeld canonical method and thus establish the non- 
perturbativc technique used in this work. 



III. THE LEWIS— RIESENFELD CANONICAL 
METHOD AND THE SINGLE PARAMETRIC 
OSCILLATOR 

In order to study the nonequilibrium dynamics de- 
scribed by time-dependent Hamiltonians, we resort to 
the Lewis-Ricscnfcld canonical method. This technique 
deals with the relation between the eigenstates of an in- 
variant operator and the solutions of the corresponding 
Schrodinger equation pTj . Using this method, it is pos- 
sible to find the exact solution for a single parametric 
quantum oscillator [ll|, [T^, [U[ • 



A. The Lewis— Riesenfeld canonical method 

An explicitly time-dependent invariant operator O 
should satisfy the condition 



1^ = 1^ 



i[H(t),6(t)]=0 



(5) 



the Schrodinger equation. In general, if \ip, t) is a solu- 
tion, this can be written in terms of the eigenstates of 
the operator 0(t), 



^c„exp(i 9 n (t))\X n ,t) , 



(7) 



where 0(t)\X n ,t) = X n \X n ,t), and the time dependent 
phase is 9 n {t) = J (A„,i|i^ - H(t)\X n ,t)dt. 



B. The single quantum parametric oscillator 



We consider the Hamiltonian 



H(t) = \f 



(8) 



where u)(t) is an arbitrary (not necessarily real), 
piecewise-continuous function of time [Tlj | . For this sys- 
tem, one can derive a set of time dependent invariants 
which are homogeneous, quadratic expressions in the dy- 
namical variables p and q. For this purpose, we introduce 
time dependent ladder operators a(t) and a'(t), that sat- 
isfy Eq. ([5]). For example, the annihilation operator is 
given by 



a(t) = -B(t)q + B(t)p 



(9) 



where the auxiliary dynamical parameter B(t) satisfy the 
classical equation of motion 



B{t)+oj 2 {t)B(t) = 
and the Wronskian condition 

B(t)B*(t)-B(t)B*(t) = 



(10) 



(11) 

The last condition ensure that the operators are ladder 
operators at all times [H, HH . In the context of Lewis- 
Riesenfeld's canonical method, we construct the Hermi- 



tian invariant operator n{t) 



it)a{t) and the Fock 



space consisting of time-dependent number states 

h(t)\n,t) = n\n,t) , (12) 

where 

(at( t ))» 



-|o,t> , 



(13) 



and the vacuum state |0, t) is anihilated by a(t). 
The dynamical behaviour of the system is fully deter- 
mined through the solution of the corresponding classical 
problem Eq. fTU]). 



We write the system's dynamics by means of the standard 
Schrodinger equation 



(6) 



A non-equilibrium system evolves towards a final state 
which, in general, differs from the initial one. Lewis and 
Riesenfeld observed that any operator 0(t) satisfying Eq. 
([5]) can be used to construct the exact quantum states of 



IV. THE LEWIS— RIESENFELD CANONICAL 
METHOD FOR COUPLED OSCILLATORS AND 
CLASSICAL PARAMETRIC RESONANCE 

In this section we find the solution of the Schrodinger 
equation for the Hamiltonian Eq. ([2]) and study the sta- 
bility properties of the auxiliary dynamical parameters, 
which arc related to the phenomenon of parametric res- 
onance in classical mechanics [3|. 
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A. Exact solution of the Schrodinger equation 

In the case of the Schrodinger equation for the Hamil- 
tonian Eq. @ , there is a transformation that allows us to 
map this Hamiltonian into a Hamiltonian of two uncou- 
pled paramctrically driven quantum harmonic oscillators. 
We define the time-independent unitary transformation 
U acting over the system's state \^,t), to obtain 



|*, t) R = U\V, t) 



(14) 



We choose the unitary operator U such that the system's 
Hamiltonian H(t) is taken to a new operator which can 
be associated to that of two uncoupled time dependent 
quantum harmonic oscillators, H R {t) = tl^H{t)U, such 
that 

H R (t) = l(f_ + (e.(t)) 2 e+pl + (e + (t)fql) , (15) 

where the energies of the two independent oscillator 
modes e T (t) are defined in terms of the other control pa- 
rameters as (ezp(i)) = ui 2 ^f2uj\j(t)\. The solution of the 
Schrodinger equation for the Hamiltonian of Eq. (|15[) is 
obtained as a linear combination of product of quantum 
states 

(t) + 8 n+ (t))\n-,t) R \n + ,t) R , 

(16) 

where \n T ,t) R are the quantum states corresponding to 
the oscillator mode with excitation energy e^(t). The 
fundamental excitations of the system are given by the 
energies e T (t). The general solution of the Schrodinger 
equation is given by 



|*,t>* 



E 



n_ ,n_j_ j 



t)i 



(17) 



In order to find the wave functions corresponding to the 
quantum states \n^,t} R , we now go back to the Lewis- 
Riesenfeld's canonical method. The first step in the con- 
struction of the solution is to find the ground state of the 
system, which must satisfy the condition 



a T \0 1 1 t) R = 



(18) 



where a=p(i) = —Bz^it)q-^ + B^(t)p^ are annihilation op- 
erators corresponding to the oscilator mode with energy 
£zp(t) that satisfy Eq. ([5]), and the auxiliary dynamical 
parameters B T (t) satisfy the Mathieu equation [l4[ 



B T (t) + (e T (t)) 2 B T (t) = , 
subject to the Wronskian condition 

B T (t)B*(t)-B^(t)B*(t)=i 



(19) 



(20) 



We next solve the differential equations subjected to an 
initial thermal equilibrium condition 0, [lfl | satisfying the 
Wronskian constraint; this requires that 



B±(t ) 



B±{t ) 



V2(e±(to))* ' 



e±(to) 



(21) 



In the coordinate representation, the operator Eq. |18p 
is equivalent to the differential equation 



-B T (t)q T + lB T (t)-^-j ^(q T ,t) = . (22) 



The solution to this equation is the wave function 



2n\B T (t)\' 



1/4 



■ B T (t) 2 



(23) 

thus obtaining the normalized ground state wave func- 
tion $§ (q_,q + ,t) = $Q (q_,t)$Q The construc- 
tion of the Fock space consisting of time-dependent num- 
ber states is direct, as any number state is obtained by 
applying the creation operators to the vacuum state 

(6 f (t))"- (at (*))"+ , 
\n-,n + ,t) R = y - -}JL K ±±=U=- \0,0,t)R. (24) 

v( n -) 1 v( n +V- 

In the coordinate representation, the number state 
is given by $%_ !n+ {q-,q+,t) = $-_{q_,t)$+ + {q + ,t), 
where the quantum states $J can ^ e written in terms 
of the Hermite polynomials H n [l3[ as 



*L(9=F.t) = 



2 -2(„ T ) \ 1/4 



2n\B T (t)\ 2 (n T y 2 J \B T (t) , 

I , (25) 



ff„(Q T (t))exp(i t(<) - 2 



, 2B=f(*) 

where we have introduced the parameter Q^(t) = 
(2\B T (t)\ 2 )~ 1 / 2 q T . In analogy to the harmonic oscilla- 
tor case, we define a characteristic system length scale in 
terms of the auxiliary dynamical parameters as 



P(t)=V2\B T (t)\ 



(26) 



The solution of the Schrodinger equation for the original 
Hamiltonian is obtained by applying the inverse of the 
unitary transformation U to the quantum state Eq. (|25|) . 
Finally, the state of the system (xi,X2\^f n -,n+,t) can be 
written as a function of the coordinates x\ and X2 as 
follows 
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0.40 



FIG. 1: Stability of the Mathieu's equation: i) shaded areas correspond to a real Floquet exponent, and the white areas to a 
complex FE; ii) the upper curve shows the real part of the FE for the solution B+(t) whereas the lower curve gives the real 
part of the FE for the solution B-(t); iii) imaginary part of the Floquet exponent for the solution B—(t) in the gap zone of the 
lower curve in ii). 



^ n -,n + (xi,X2,t) = exp(i 6 n _ (t)) exp(i 9 n (t))$ n _ 



r 



V2 yft'J " + lV2 + V2' 



(27) 



B. Mathieu's equation and parametric resonance 



Mathieu's equation [3] is given by 

f(t) + (a - 2fccos(2i))/(i) = 



(28) 



the behaviour of its solution depends strongly on the pa- 
rameters a and b, the driving can stabilize or destabilize 
the undriven oscillation, this means that the solution can 
be bounded or increasing with time respectively. In gen- 
eral, parametric resonance can occur if the parameters of 
a classical dynamical system vary periodically with time. 
Stable fixed points of the flow in phase space become 
unstable for specific values of certain parameters 0] . In- 
terestingly, the set of differential equations (|19p has the 
general form of Mathieu's equation [l4| 

B T (t) + [(lj 2 =f 2ug) T (2wA 5 ) cos ilt]B T (t) = . (29) 

If we compare this with the canonical form of the Math- 
ieu's equation of Eq. (|28|) , it is clear that the parameters 
a and b are functions of the static coupling g, the fre- 
quency w, and the external driving frequency f2: 



b T (g,u;,n) 



A(oj 2 T2ujg)/il 2 . 
±4u;Ag/fl 2 . 



(30) 
(31) 



In order to explore the relevant characteristics of the so- 
lutions of Eq. (|2"9"|) , we apply the Floquet theorem for 



second order differential equations with time periodic co- 
efficients [lH . For this, the solutions of Eq. (|29p have the 
general form 



B T (t) = exp(iF T ^ T (*) 



(32) 



where (f> T (t+T) = <j> T (t) and F T (a T (g,uj,n),b T (g,uj,n)) 
is the Floquet exponent (FE), which depends on the 
shape of the driving. The stability zones of the Math- 
ieu's equation are well known: for driving functions for 
which Fzf: is complex, the solution becomes unstable, and 
in the stable regime, F T is a real number [3, Ell- ^ n 
Fig. Q] (i) we show the stability zones of the Mathieu's 
equation in the a — b plane: within the shaded areas the 
FE is real, whereas within the white areas this becomes 
complex. Once we have established the parameters Ag, 
oj, and f2, the map C T (9) = ( a=F (5? w ? b T (g,uj, f2)) de- 
scribes a straight line in the a — b plane, as can be ob- 
served in Fig. [T] (i). The real part 3?(F T (C T (g))) of the 
FE for the auxiliary dynamical parameters B^(t) as a 
function of the static coupling, exhibits band-like regions 
that correspond to a complex FE of the solution, which 
means that this solution is unstable within the "gaps" , 
where the Floquet exponent is imaginary. In order to 
study such features, we first investigate the particular 
case Ag = O.lg, u = 1, and fi = 1. Fig. Q] (ii) depict 
the real part of the FE for the solution B T (t) and Fig.Q] 
(iii) the imaginary part of FE for the solution B~(t). In 
these figures, we find that for this particular case only ex- 
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FIG. 2: Phase space trajectories of the auxiliary dynamical parameters B (t) for A<? = O.lg, and a value of the static coupling 
in a) the unstable region, g = 0.38, b) the stable region, g = 0.4, and c) the stable region, g — 0.46. 



ist one band-like region for the solution B~ (£). In Fig. [21 
we plot the phase space representation of the trajectories 
of the auxiliary dynamical parameters B T (t). Figure [2] 
(a) shows that when the static coupling belongs to the 
instability zone (g = 0.38) of B_(t), the solution B + {t) 
is bounded, while the solution i?_ is unbounded in the 
phase space. In contrast, in Figs. [2^6) and [2] (c), it is 
clear that for values of g that belong to a mutual stability 
zone (g = 0.4, and g = 0.46) the solutions are bounded 
in the phase space. 

We can find an exact solution to the Schrodingcr equation 
for the general problem described by the Hamiltonian Eq. 
([Tl) through the Lewis and Riescnfcld canonical method 
[ill . H2I [l3l | , the general solution is given in Appendix A. 



V. QUANTUM DYNAMICS AND 
ENTANGLEMENT QUANTIFIERS 

In this section we establish a connection between the 
non-equilibrium entanglement generated by the time- 
dependent coupling and the behaviour of the auxiliary 
dynamical parameters. Furthermore, in order to study 
the correlations in the ground state and their relation 
with the localisation of the total system's wave function, 
we derive analytical expressions for the reduced density 
operator of one oscillator and the corresponding Wigner 
function. 



A. The reduced density operator 

When studying a composed quantum system whose dy- 
namics is unitary, it is often more interesting to study a 
subsystem of the whole system, which is described by 
its reduced system dynamics. In contrast to the dy- 
namics of the whole system, the dynamics of the sub- 



system is not unitary. It is possible, however, to make 
such a description of the subsystem through the reduced 
density operator. We describe the dynamics of the to- 
tal system in terms of the density operator, where the 
pure state of the universe is represented by the operator 
Pa{t) = |^o_, o+i t)($fo-, o+) t\- In the coordinate repre- 
sentation, the density matrix takes the form 

p G (x' 1 ,x 2 ;x 1 ,x 2 ,t) = % Q+ (x 1 ,x 2 ,t)^ ^ t0+ (x' 1 ,x 2 ,t) . 

(33) 

At this stage, we develop a bipartite decomposition of the 
universe, which lets us study the reduced dynamics of one 
oscillator through the reduced density matrix (RDM) . In 
so doing, we calculate the partial trace over the physical 
field coordinate x\ , 

/+00 
*0_ ,0+ { x l 1 X 2 , *)*<)_ ,o+ {xi ,x' 2 ,t) dxi . 
-00 

(34) 

The ground state of the universe ^o_.o + (xi,x 2 ,t) given 
in Eq. (|2"T)) is Gaussian, which facilitates the integration 
ofEq. (f3~4"]) . After some algebraic calculations, we obtain 
the following result for the reduced density matrix 

p^\x' 2 ,x 2 ,t) = Acxp(^(a)((x' 2 ) 2 + x 2 2 ) + l3x 2 x' 2 )x 
exp(i%(a)((x' 2 ) 2 -x 2 2 ))), (35) 
where the parameters a, j3, and A are defined as 



'[(r-£ + )' 



(n*s 2 + (£ + )*c 2 

2 4(5R(£-)c 2 + 9?(£+)s 2 )' 

cV(r-e + )*(r-e + ) 

2(»(£-)c 2 + 9fi(£+)s 2 ) ' 



A = 



1/2 



(36) 



wm~)c 2 +m + > 2 ) 

where c = s = l/\/2 and 5i(£ T ) is the real part of the 
function £ T (t) = — i a (t) • ^ n tne particular case j(t) = 



6 




g, this function is independent of time and becomes 



(37) 



and corresponds to the result previously reported in the 
literature for the atom's reduced density matrix for the 
single mode Dicke model in thermal equilibrium [^, [To| . 
The reduced density operator Eq. (|35j) has the same 
structure as the density operator of an ensemble of time- 
dependent oscillators in the coordinate representation 



elepe 
0- 



B. The linear entropy 



A density operator p describes a pure state if and only 
if it satisfies p 2 — p, i.e., the density operator is idem- 
potent. We consider a pure state of a composed system 
(the universe), say system AB, represented by the den- 
sity operator pab- The linear entropy for the reduced 
density operator pa = trsipAB) of subsystem A is de- 
fined by La = 1 — trv^p^) [lj. This gives a measure 
of purity of the reduced density operator pA (one part 
of the total system, in a bipartite decomposition of the 
universe). If the pure state of the universe is separable, 
the reduced density operator of one part of the system 
represents a pure state and, as a result, its linear entropy 
must be zero. Similarly, if the pure state of the universe 
is a maximally entangled state [13], the linear entropy 
equals 1/2. In the context of our model, it is possible to 
use the linear entropy as a measure of bipartite entan- 
glement for the oscillators system. The linear entropy 
for the oscillator's time dependent reduced density oper- 



ator, Eq. (|35[) . is given by L(t) = 1 — tr 



By using the density matrix representation of Eq. 
we explicitly calculate the linear entropy as 



L(t) = 1- 



— OO J — OG 



/°g ?) ( x 2 ' x 2 , t)p ( a' {%2 ,x' 2 ,t) dx 2 dx' 2 , 



„(*) 



which gives, after some algebra, the result 



L(t,g) = 1 



ttA 2 



^4(5R(a)) 2 - P 2 



(38) 



In Fig. [3]we plot this quantity as a function of time: it is 
clear that in the case of a time periodic atom-field cou- 
pling there exist a direct correspondence between the sta- 
bility of the auxiliary dynamical parameters and the en- 
tanglement dynamics, which is reflected in the graph[3la) 
for a value of the static coupling in the unstable zone, 
it is observed that the linear entropy oscillates, before 
reaching the stationary state with diverging linear en- 
tropy L = 1. We consider the linear entropy in Figs. [3] 
(b), and[3](c) for values of the static coupling in the sta- 
ble zone. Our choice of the initial conditions for Eq. (TT5|) 
implies that the system starts in thermal equilibrium, 
but this is described by a non-separable quantum state. 
The dynamics exhibits a behaviour whereby the system 
experiences a disentanglement process, before reaching 
a maximum value of entanglement, and then successive 
collapses and revivals. The frequency of the oscillations 
in the "stable zones" depends strongly on the value of 
the static coupling g. 
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C. The Wigner function 

The Wigner function is a particular representation of 
the density operator for a pure or a mixed state of a 
quantum system. In principle, the density operator is 
more fundamental than its Wigner representation, how- 
ever, the Wigner function is often a useful tool in deco- 
herence studies for investigating possible correlations be- 
tween position and momentum [f| [lU, [3 . The Wigner 
representation of a given density operator p is given by 



1 f°° 

W{q,p) = — / exp (ips) (q- 



s/2\p\q + s/2) da, (39) 



where we have considered h = 1. The Wigner function 
can take on negative as positive values. This means that 



J 



it cannot be strictly interpreted as a phase space proba- 
bility density. In the special case when the density opera- 
tor represents a pure Gaussian state, the Wigner function 
is Gaussian and hence is positive definite. The reduced 
Wigner function represents the partial trace of a density 
operator over a subsystem and so contains all informa- 
tion about a given subsystem. Here we calculate the 
Wigner representation of the oscillator's reduced density 
operator given by Eq. (|35[) . Such Wigner function reads 



1 

2^ 



e lpa p { cXq- s/2,q + s/2,t) ds 



By performing the Gaussian integral we obtain 



(40) 



W (R) (q,p,t) = A 1 



1 



7r(23?(a) + (3) 



exp 



43(a)gp 




exp 


_29?(a) +/?_ 



g 2 (4|q| 2 -/? 2 )+p 2 
2SR(a) + (3 



(41) 



r 



where the parameters a, j3, and A are defined as in Eq. 
(|36|) . If we specialize to the case j(t) = g, the Wigner 
function becomes time-independent, 



W {R) (q,p) 



?r(e_ +e+) 



exp 



P 



(42) 

where arc defined in Eq. (|37[) . This function exhibits 
some interesting features: for values of the coupling near 
to the critical coupling g c — ^, the Wigner function be- 
comes stretched along the position axis with a consequent 
contraction along the momentum axis in the phase space; 
this results in the Wigner function being delocalized in 
phase space. This behaviour is due to the existing entan- 
glement, because the oscillator's state becomes a mixed 
state. Such a behaviour is corroborated by a divergence 
of the von Neumann entropy near the critical coupling, 
signaling a "maximal mixing" scenario, as reported in 
Ref. [T(| in the context of the time-independent Dickc 
model. Figures 2] and [5] show the evolution of the Wigner 
function in the case of a periodic atom-field coupling for 
values of the static coupling that belong to unstable and 
stable zones, respectively. In our dynamical coupling sce- 
nario (Ag 7^ 0), there is a direct relationship between 
the atom-field entanglement and the localisation of the 
Wigner function in phase space: the Wigner function 
evolves in a delocalised region of the phase space for val- 
ues of the static coupling in the unstable zone; as shown 
in Fig. [31 the function evolves in such a way that it is 
stretched along a dynamical rotating axis with a conse- 
quent contraction in the perpendicular axis, in contrast 
to the time- independent behaviour. The behaviour ex- 
hibited by the Wigner function for a value of the static 



coupling in the stable zone is quite different: the Wigner 
function also becomes stretched along a dynamical axis 
but it is now localised as the system evolves, as shown in 
Fig. El 



VI. DISCUSSION 

From our analytical calculations, we have found that 
in contrast to the time independent case, the linear en- 
tropy exhibits a time dependence that is not determined 
by the global phase of universe's ground state wave func- 
tion, Eq. I|27p. but by the auxiliary dynamical param- 
eters Bzf(t). Numerical results suggest that when the 
set of auxiliary dynamical parameters are under para- 
metric resonance, the trace of (p^\t)^ decays expo- 
nentially with time with periodic modulations. A care- 
ful exploration of the instability zone, shows that the 
imaginary part of the Floquet exponent grows from zero 
within a certain interval and then falls to zero as shown 
in Fig. [T]Jm). This fact involves a transition between os- 
cillatory and diverging linear entropy and is related with 
the upper quantum Lyapunov exponent pol . |2lj . The 
characteristic length of Eq. (|2l))) is bounded in the stable 
zones and unbounded in the unstable zones. The prod- 
uct of the two characteristic length scales in the system is 
proportional to the Gaussian normalisation factor of the 
ground state, a result that tells us about the relative vol- 
ume in coordinate space that the wave packet occupies. 
The numerical simulations of the ground state probabil- 
ity density show that for values of the static coupling 
that belong to the common stability zones, the proba- 
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FIG. 4: Time evolution of the Wigner function in p — q space for the "unstable zone" g = 0.38, and Ag = O.lg; a) Q,ti = 0, b) 
nt 2 = 32, and c) nt a = 50. 




FIG. 5: Time evolution of the Wigner function in p — g-space for the "stable zone" g = 0.4, and Ag = 0.1; a) Qt\ — 0, b) 
Q.t 2 = 32, and c) tlt s = 50. 



bility is localized in the x\ — space and exhibits an 
oscillatory behaviour. In contrast, for values of the static 
coupling that belong to the unstable zones, the probabil- 
ity density has also an oscillatory behaviour but with a 
different ingredient: when the system evolves, the den- 
sity is systematically stretched in a fixed direction, with 
a consequent dilation in the perpendicular direction, this 



behavior results in a massive delocalisation of the wave 
function of the total system. This behaviour is very sug- 
gestive, because the Wigner function for one oscillator 
in our case exhibits a similar pattern, but with a cor- 
responding change in the stretching direction in phase 
space as time evolves, the Wigner function of one oscil- 
lator ensemble results in a delocalisation in phase space, 
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a fact that we associate with the "maximal mixing" sce- 
nario. 

Recently, Blume-Kohout and Zurek analyzed a 

model in which they considered a harmonic oscillator 
coupled to an unstable environment, which consisted of 
an inverted harmonic oscillator. They demonstrated that 
this unstable environment could produce decoherence in 
the system of interest more effectively than a bath of har- 
monic oscillators. They used the von Neumann entropy 
as a measure of bipartite entanglement, and surprisingly, 
entropy increases linearly over time with periodic modu- 
lations. In contrast to this result, when the environment 
is a stable oscillator, the von Neumann entropy oscil- 
lates, but does not increase over time. We can think of 
our system from another perspective: assuming that the 
universe is composed of two coupled oscillators, we con- 
sider that one oscillator is the system of interest, and the 
other oscillator is the environment. In contrast to the 
model 22], in our case, both systems are stable, how- 
ever, control over the dynamics is established through 
the time-dependent coupling between the oscillators. Nu- 
merical results suggest that the von Neumann entropy for 
the system's reduced density operator exhibits the same 
behavior described above. Comparing our results with 
those found by Blume-Kohout and Zurek, we find that 
the effect of external control over the system is to pro- 
duce instability in the universe, which is related to the 
phenomenon of parametric resonance. This instability is 
reflected in a loss of coherence in the system of interest. 



VII. CONCLUSIONS 



We have obtained exact results for the oscillator's re- 
duced density operator and the linear entropy. The 
results allow us to study the entanglement dynamics 
through the calculation of the auxiliary dynamical pa- 
rameters. Our results for the Wigner function are to be 
contrasted with the time independent case, where the en- 
tanglement does not exhibit dynamics because the tem- 
poral dependence of the ground state of the universe 
is given by a global phase of the wave function. We 
have found that near the critical coupling, the oscillator's 
Wigner function in the time-independent case results in 
a dclocalisation in phase space, a fact that we associate 
with the "maximal mixing" scenario reported in (Io| via 
the divergence of the von Neumann entropy. We have 
given a detailed prescription for the new entanglement 
features associated to the system in terms of the linear 
entropy and its relation to the Wigner function. The aux- 
iliary dynamical parameters and its stability properties 
play a crucial role when determining the entanglement 
properties of the system and the behaviour of the uni- 
verse ground state wave packet. 
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APPENDIX A: SOLUTION OF THE 
SCHRODINGER EQUATION FOR A 
TIME-DEPEDENT QUADRATIC 
HAMILTONIAN 

We use the Lewis-Riescnfcld canonical method to find 
the exact solution of the Schrodinger equation for a sys- 
tem described by the Hamiltonian Eq. |T]). With the 
purpose of building the invariant operator, we define the 
time-dependent ladder operators 

2 

= J2(A ik (t)x k + B ik (t)p k ) , (Al) 
fc=i 

and 

2 

al(t)=Y / (A* k (t)x k + B* k (t)p k ), (A2) 
fc=i 

which satisfy the standard bosonic commutation rela- 
tions 



[a m (t),a n (t)} = [at i (t),6+(t)]=0 1 



(A3) 



and 



[a m (t),al l (t)]=6 mn . (A4) 

The operators a,(t) and aj(f) must satisfy the equation 
Eq. ([5]), a condition which implies that the auxiliary 
dynamical parameters Ai k and Bi k satisfy the following 
differential equations 

Bu(t) = ^StBh -A*i(t)*l(*)-B»-i(*) 

-^(t) 7 (t)B i2 (t) , (A5) 

B i2 (t) = ^B a -p 2 {t)v 2 {t)B l2 {t) 

-/i 2 (t)7(t)fl i i(t) , (A6) 

and 

B ik = -fM k (t)A ik (t) . (A7) 

As a consequence of these requirements, the auxiliary 
dynamical parameters B ik (t) should satisfy 



B mk (t)B nk (t) 


B nk {t)B mk (t) 


Mfc(*) 


Hk{t) 


B* mk (t)K k (t) 


K k (t)B* mk (t) 




M2(*) 


B m k(t)B* k (t) 


B mk (t)B* nk (t) 



p k (t) 



,(A8) 



= ,(A9) 



»<WAio) 
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Using Eqs. (|A5[) and (|A6[) . one can easily show that Eqs . 
(HHD , ([AT)]) , and (jATU)) must be time independent [llEl. 
Interestingly, the auxiliary dynamical parameters Bo- (t) 
are solutions of the classical equations of motion for a 
non-conservative Hamiltonian system 0, [l3[ , as is clearly 
seen in Eqs. (|A5|) and (|A6j) . This fact is related to the 
Ehrcnfest theorem, because for a particle in a parabolic 
potential well, the motion of the center of the wave 
packet rigorously obeys the laws of classical mechanics 
[23[ ■ The previous procedure allows us to define the Her- 
mitian invariant operator 0(t) = hi(t) + 77.2 (i), where 



hk(t) = a k (t)ak{t) are time-dependent number opera- 
tors. We construct the Fock space consisting of time- 
dependent number states |ni,ri2,t), where the vacuum 
state |0, Q,t) is the only one that is annihilated by &i(t) 
and 02 (i). The number state is obtained by applying 
a\(t) and a 2 (t) to the vacuum state 

\n l ,n 2 ,t) = ^Mm\ W ). (All) 
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